#include "../Random.hpp"
#include <armadillo>
#include <boost/math/distributions/normal.hpp>
#include "../utils.hpp"


#define _Mx 100

template<class Approx,class APriori,class Likelihoodms>
class Probit_MS2 : public Density::GeomBridge
{
		friend class Density::GeomBridge;
	public:
		Probit_MS2(double *Y, mat X, Distribution::Distribution *S,APriori *P) : Density::GeomBridge(Y, X, S){ 
			_s=P->Get_s();
			_d=new boost::math::normal_distribution<>(0,1);
			Set_p(X.n_cols-1);
			_is=inv(_s);
			_m=_Mx;
			fL=new Approx;
			_Px=P;
		}
		~Probit_MS2(){
			delete fL;
		}
		inline double Likelihood(mat gamma)
		{
			int p=Get_p();
			mat X=Get_X();
			mat gammat(p+1,1);
			mat h=sum(gamma,0);
			double res=0;
			int nb=h(0,0)+1;//nb column (use for filter)
			int n=X.n_rows;	
			mat X2(n,nb);
//			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(p+1);i++)
			{
					
				if(i!=0){
					gammat(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=X(span::all,i);
						kk++;
					}

				}else{
					gammat(i,0)=1;
					X2(span::all,0)=X(span::all,0);
				}
			}
		
			mat Xt=X*diagmat(gammat);
			mat m(nb,1);
			m.zeros();
			double *Y=Get_Y();
		//	cout << X2;	
			(*fL)(X2,Y,m,_s(0,0));			
			mat Sig=(*fL).Get_Sig();	
			mat Mu=(*fL).Get_Mu();
			///////////////////////
			mat bar(nb,nb);
			bar.fill(0.1);
			mat S=diagmat(bar);
			bar.fill(10);
			_s=diagmat(bar);
			G=new Distribution::Gaussian(nb,m,S);
			F=new Distribution::Gaussian(nb,Mu,Sig);
			DP=new Likelihoodms(Y,X2,F,_s);
			//////////////////////
			K=new Kernel_MH<Density::GeomBridge,Distribution::Gaussian>(DP,G);
			//Distribution::Gaussian G(Xt.n_cols,m,s);
			R=new Resample::Multinomial;
			Particle<Kernel_MH<Density::GeomBridge,Distribution::Gaussian>,Resample::Multinomial,Density::GeomBridge> P(K,R,DP,_Mx,F,0.5);
			P.Filter();
			//cout << "test" << nb-prop.n_rows;
			double sum=P.Get_Z();	
			//cout << res << "\n";
			delete F;
			delete DP;
			delete G;
			delete R;
			delete K; 

			return sum;

		}

		double Prior(mat theta){
			return -0.693;
		}
		mat  GradLik(mat theta)
		{
			cout << "Grad ProbitMS deprecated";
			return theta;
		}

	private:
		
		Kernel_MH<Density::GeomBridge,Distribution::Gaussian> *K;
		Distribution::Gaussian *G;
		Distribution::Gaussian *F;
		Distribution::Distribution *_Px;
		Resample::Multinomial *R;
		Likelihoods::Probit3 *DP;
		Laplace *fL;
		boost::math::normal_distribution<> *_d;
		mat _s;
		int _m;
		mat _is;
	
};
